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Two-dimensional gyrokinetics is a simple paradigm for the study of kinetic magnetised 
plasma turbulence. In this paper, wc present a comprehensive theoretical framework for 
this turbulence. We study both the inverse and direct cascades (the 'dual cascade'), 
driven by a homogeneous and isotropic random forcing. The key characteristic length of 
gyrokinetics, the Larmor radius, divides scales into two physically distinct ranges. For 
scales larger than the Larmor radius, we derive the familiar Charney-Hasegawa-Mima 
(CHM) equation from the gyrokinetic system, and explain its relationship to gyrokinetics. 
At scales smaller than the Larmor radius, a dual cascade occurs in phase space (two 
dimensions in position space plus one dimension in velocity space) via a nonlinear phase- 
mixing process. We show that at these sub-Larmor scales, the turbulence is self-similar 
and exhibits power law spectra in position and velocity space. We propose a Hankel- 
transform formalism to characterise velocity-space spectra. We derive the exact relations 
for third-order structure functions, analogous to Kolmogorov's four-fifths and Yaglom's 
four-thirds laws and valid at both long and short wavelengths. We show how the general 
gyrokinetic invariants are related to the particular invariants that control the dual cascade 
in the long- and short-wavelength limits. We describe the full range of cascades from the 
fluid to the fully kinetic range. 



1. Introduction 

1.1. Background 

A fluid is conventionally described by macroscopic state variables such as bulk flow ve- 
locity, density, pressure, etc., which vary over three-dimensional space. This description 
is appropriate when collisions between the constituent particles establish local thermody- 
namic equilibrium (Maxwellian velocity distribution) more rapidly than any dynamical 
processes can disturb this equilibrium. When this condition is not met, a kinetic descrip- 
tion is needed to capture the evolution of a distribution function in six-dimensional phase 
space (positions and velocities). 

Weakly collisional plasmas require a kinetic description to capture a wealth of dynam- 
ical phenomena that are wrapped up in the velocity distribution of particles. From the 
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perspective of a turbulence theorist, the foUowing questions are particularly interesting: 
What are the dynamical invariants? How do they travel through the phase space? How 
do the velocity dimensions participate in the fully developed state of the turbulence? 
Can one define an "inertial range" for the phase space in kinetic turbulence? How can 
one, in general, adapt the phenomenology of fluid turbulence — an energy cascade by 
local interactions, universal self-similar scaling, dissipation scale, and so forth? 

In this paper we address these questions for a simple kinetic system: the two-dimensional 
electrostatic^ gyrokinetic system driven by a statistically homogeneous and isotropic 
source. The gyrokinetic system of equations ( Taylor fc Hastie||1968 [Rutherford fc Frie 



man'196 8| | Catto fc Tsang|1977|[Antonsen fc Lane|1980||Frieman fc Chen|1982||Brizard fc 



Hahm 2007 ) is used to describe magnetised plasma dynamics on time scales much larger 
than the ion Larmor periodj^ In this limit, the presence of a strong, constant magnetic 
"guide" field causes the charged particles to rapidly trace nearly circular Larmor orbits so 
that the electrostatic fiuctuations affect the motion adiabatically. This causes slow drift 
motion perpendicular to the guide field. Thus, what we are considering is essentially a 
kinetic theory of charged "rings." 

We assume that the distribution function for the particles is a Maxwellian distribu- 
tion plus a perturbation. As a result, the intrinsic velocity scale in this problem is the 
thermal velocity, Wth- Because of the magnetic guide field, this velocity scale gives rise 
to an associated intrinsic spatial scale which is p, the thermal Larmor radius p — wth/^^, 
where f2 is the Larmor frequency. Physically, the kinetic phenomenon on which we will 
be focusing arises as follows. For electrostatic fluctuations on scales much larger than the 
Larmor radius, all particles move together with the same E x B drift and a fluid descrip- 
tion is correct. At the Larmor-radius scale and smaller, particles of distinct velocities 
have different effective E x B drifts (as they sample different regions of the electric field 
during their rapid Larmor motion) and a kinetic description is required. At such scales, 
gyrokinetic turbulence exhibits a distinctive kinetic behaviour in phase space due to the 
presence of so-called "nonlinear phase- mixing" ( Borland fc Hammett|1993 Schekochihin 



\et aL 2008 Tatsuno et a/. ^ 2009) . Fluctuations in the distribution function nonlinearly 



cascade to create fine-scale structure in velocity space in addition to the two real-space 
dimensions — thus, we may think of velocity space as an additional dimension to be 
treated on equal footing with position spacej^ 

To address the problem of gyrokinetic turbulence, we borrow traditional methods from 
fiuid turbulence. This is possible because the gyrokinetic system bears some resemblance 
to equations of incompressible fluid turbulence. The nonlinearity is convective, with a 
divergence- free velocity (the E x B drift velocity due to fluctuating electrostatic flelds), 
while dissipation occurs via a collision operator that acts by second-order derivatives (see 
Catto fc Tsang|[T977 Abel et aL]|2008 ) in phase space and, therefore, acts on fine scales 



in analogy to viscosity in fluid turbulence. In the long-wavelength cold-ion limit (see sec- 
tion l4|, the electrostatic gyrokinetic system can be reduced to the Charney-Hasegawa- 



^ The electrostatic limit in this context means that there are no magnetic fluctuations, al- 
though there is a strong mean magnetic field, constant in both space and time. 

^ Gyrokinetics has been developed by the magnetic fusion community to describe plasma 
turbulence that causes the transport of heat and particles in fusion de vices. A review o f the 
theoretical framework of turbulence in magnet ised plasmas is given by 'Krommes' (' 2002]). As 
demon strated in the series of recent works by iHowes et al. ( 2006, 2008); Schekochihin et al. 



(20091, gyrokinetics is also appropriate for a wide range of astrophysical plasmas. 

Note that the nonlinear phase- mixing studied in this work depends critically on the assump- 
tion that the plasma is magnetised. One should not necessarily expect the type of phase-space 
cascade considered here to appear in all types of kinetic turbulence. 
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Mima (CHM) equation fCharney' |1971[ |Hasegawa fc Mima|[l978[ ) or, with additional 
assumptions, to the inviscid vorticity equation that describes two-dimensional Euler tur- 
bulence (Kraichnan 1967 Batchelor 1969[ )P] Thus, from a fluid turbulence perspective, 
two-dimensional electrostatic gyrokinetics can be seen as a simple kinetic extension of 
extensively-studied fluid equations. 

1.2. Outline 

The progression of the paper is as follows. We begin in section [2] by discussing the 
gyrokinetic system, its applicability and the specific assumptions and approximations we 
will be using. In section[3| we introduce and discuss the dynamical coUisionless invariants 
of the gyrokinetic system. It is shown that there are two such invariants. One invariant 
is related to the perturbed free energy (or perturbed entropy) of the system. We argue 
that this invariant must cascade forward to fine scales in phase space, where the collision 
operator ultimately dissipates it. The other invariant is special to the two-dimensional 
electrostatic system and is not conserved in three-dimensional gyrokinetics. We argue 
that this "electrostatic" invariant will cascade inversely to larger scales — as does energy 
in two-dimensional fluid turbulence. 

In section |4] we study the long wavelength regime (fcp <^ 1) and explore the relation- 
ship between gyrokinetic and CHM/Euler turbulence. We derive the inviscid CHM/Euler 
equations as asymptotic limits of the gyrokinetic system. We also discuss the two in- 
variants of the CHM/Euler equation and how they relate to the gyrokinetic collision- 
less invariants. We find that our electrostatic invariant in this limit becomes the in- 
versely cascading CHM/Euler invariant ('energy'). The forward cascading quantity from 
CHM/Euler turbulence ('enstrophy') is found to be only one part of the gyrokinetic quan- 
tity which cascades forward; the other essentially kinetic part corresponds to a passive 
scalar field. 

In section [Sj we derive several scaling results by phenomenological arguments. This 
section gives a simplified physical perspective, which may serve as a guide for the analysis 
in the remaining sections of the paper. We begin with a standard treatment of the long- 
wavelength fluid regimes. We then shift attention to the short-wavelength {kp 3> 1) 
nonlinear phase-mixing regime, where the turbulence has a strongly kinetic character — 
this regime is the chief focus of the remainder of the paper. We describe the physical 
mechanism of nonlinear phase-mixing and present a phenomenological analysis of the 



turbulence cascade based on the work by Schekochihin et al. (2008 2009). In section 5.4 



we tie together the results by presenting two scenarios - injection at a long-wavelength 
scale and injection at a short-wavelength scale - and describing the dual cascade through 
various scale ranges. 

In section [6] we begin a more formal analysis of the gyrokinetic system. Wc first intro- 
duce the statistical tools and notation that will be needed, followed by a discussion of 
symmetries. Symmetries play an important role in turbulence theory, and in subsequent 
arguments we appeal to the principle of restored symmetry in the fully developed state. 
In sections |6.3| and |6.4[ we derive the third-order statistical laws following from the gy- 



rokinetic equation, in the style of Kolmogorov ( 1941 ) and Yaglom ( 1949^ . These are exact 
results, valid uniformly for all wavenumbers, i.e. both in the long- and short- wavelength 
limits. 

In section [8] we derive scaling laws for the phase-space spectra. Scales in real space are 
treated in the conventional way by using a Fourier transform while for scales in velocity 



* The relationship between plas ma dynamics and two-dimensional fluid turbulence was 
demonstrated first by .Taylor & McNamara (19711. 
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8.1. Self-similarity 



5.4. Multi-scale cascade scenarios 



.2. Spectral scaling in phase-space 



Figure 1. Structure of the paper. Unidirectional arrows indicate dependencies and flow of 
discussion. Bidirectional arrows indicate cross-referencing. 



space, we use the zeroth-order Hankel transform, which is inspired by the mathematical 
structure of gyrokinetics. Using this formalism, we derive approximate spectral scaling 
laws, both for the forward and inverse cascades, in the limit of spatial scales much smaller 
than the Larmor radius and velocity scales much smaller than the thermal velocity. 



1.3. Structure of paper 

To help the reader navigate the paper, we have provided a diagram in figure [T] that 
describes its overall structure. Although the paper is arranged to be read from start to 
finish, the reader may choose to narrow his/her attention by following only one of the 
three suggested paths in the diagram. These paths are sufficiently independent to be read 
separately, but are strongly related. The left path in the diagram focuses on the familiar 
CHM and Euler fluid limits in the context of gyrokinetic theory. In the center path, 
many of the short-wavelength kinetic results are derived by simple phcnomenological 
arguments. The right path focuses on derivation of the results using exact methods and 
a novel spectral treatment of velocity space. There is a significant amount of material 
in sections |6) [7] and [8] which is not covered in the phcnomenological section [5] but the 
derivations are more technical and the central path should provide a good introduction 
to the important ideas. 



2. Equations 

What follows is a brief discussion of the gyrokinetic equations, their normalisation and 
applicability. A central point that emerges from this discussion is that the system studied 
in this paper will have wide-ranging applicability, subject to specific interpretation of the 
fluctuating fields and normalisation. Readers wishing to skip these details are advised to 
take equal ions |2.5| and |2.6| with definitions |2 . 3[ [2^ and [277| as the equations we will study 
in this paper. 

In kinetic plasma theory, each species s (electrons or ions) of a plasma is described by 
a distribution function /, which is a function of space, time and velocity. As is customary 
with so-called 5f gyrokinetics, the full distribution function is split into an equilibrium 
part and fluctuating parts: 
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f = Fo + e[~ipFo + h{t,IL,v^,v\\)] , (2.1) 

where if is the electrostatic potential and e = p/L is the fundamental gyrokinetic ex- 
pansion parameter, the ratio of the Larmor radius to the system scale. The equilibrium 
distribution Fq — (27r)~^/^e~^ is a Maxwellian distribution in velocity space. The 
second term in equation |2.1| is the so-called Boltzmann part of the distribution. The 
final term is the gyro-centre distribution function, h, which depends on coordinates that 
are defined relative to the mean magnetic field, Bq = zBq: the perpendicular velocity 
v± (magnitude of the velocity perpendicular to Bq), the parallel velocity v\\ (velocity 
parallel to Bq) and the gyro-centre position R. The particle position and gyro-centre 
position are related by the transformation r = R -I- p where the Larmor radius vector is 
p{'d) = z X V = wj, (y cos ■(? — X sin 1?) and ^ is the gyro-angle. Velocities are normalised to 
the thermal speed Uth, while spatial scales are normalised to the Larmor radius p = vn^/il. 
The normalisations are given in full in the appendix |Aj 

It can be shown by rigourous asymptotic expansion {e.g., see Howes et aZ.||2006 ) that 



the gyro-centre distribution function h satisfies the nonlinear gyrokinetic equation. The 
electrostatic gyrokinetic equation for a homogeneous equilibrium (i.e., assuming no vari- 
ation in the background magnetic field, density or temperature) is 

f +^„g+v..V.-(C[/.]),= ^F„ + ^, (2.2) 

where V = d/dR. The term {C[h])-^ represents the gyrokinetic collision operator, and T 
is a general kinetic source term (i.e. external drive), which will be left unspecified at this 
stage. The gyro-average is an average over the gyro-angle, with gyro-centre position 
R held fixed: for an arbitrary function of position A{r), 

(^W>R=^^'"dM(R + p(^)). (2.3) 

Note that, even if A{r) has no velocity dependence, its gyro-average (A(r))j^ does. The 
gyro-averaged E x B velocity, a function of gyro-centre position R and velocity vj_, is 
defined 



Vb(R,i;^) -z X V((^)j^. (2.4) 
In our two-dimensional theory, we neglect the parallel ion inertia term v^^dh/dz in 



equation 



2.2 



see section 



9.4 



for a discussion of this). The U|| dependence of h is therefore 
no longer important (except for its appearance in the coUisional operatoij^. Thus, it will 
be concealed when possible in the remainder of this paper. It should be assumed that 
wherever velocity integration is present, integration over t;|| is implied. Henceforth, we 
use V as shorthand for the perpendicular velocity v±. 

We can rewrite the gyrokinetic equation |2 . 2| compactly in terms of a new distribution 
function g = h — {(p)-^Fq. This is two-dimensional gyrokinetic equation that we will use 
in the rest of this paper: 



When written explicitly, t he full coUisional o perator is the gyro-averaged linearised Landau 
operator (see the appendix of [Schekochihin et al. 2009), an integro-differential operator with 
non-trivial depen dence on space. However, a simplified operator, such as the one given by 



Abel et al. 



(20081, is a good example to keep in mind for the discussions of this paper. 
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^^+^^.Vg^{C[h])^ + F. (2.5) 

The full plasma dynamics are described by the gyrokinetic equations for each species, 
with the additional constraint that charge neutrality must be maintained ( "quasi-neutrality" ] 
We make a common simplification in solving the gyrokinetic equation for only one of the 
species. This is justified by the disparity between the ion and electron Larmor radii, which 
leads to the Boltzmann (or 'adiabatic') response for the other species (see appendix [b] 
for a discussion). Under these assumptions, quasi-neutrality takes the form 

/"OO 

2tt vdv{g)^ = {l + T)if-TQ^. (2.6) 







The constant r depends on the choice of kinetic species (or scale of interest) and response 
model for the non- kinetic species (see appendix |b]). Note that the velocity integration at 
constant real-space coordinate r requires a change of variables (from gyro-center coordi- 
nate R). This gives rise to the angle average {.)^, which is an average over gyro-angle "d 
with real-space coordinate r held fixed: 

{A{K))^ = ^ j'J dm^^pm- (2.7) 

The operator Fq is defined 

poo 

T^^= / vdve-^''^{^)^)^. (2.8) 
This operator is multiplicative in Fourier space — it is simply 



fo(fc)=/o(fc2)e-\ (2.9) 
where /q is the zeroth-order modified Bessel function. Note that the gyro-average in 



Fourier space reduces to multiplication by the Bessel function Jo(fcw). So, equation 2.4 is 



= iz X kJo(few)<^(k), (2.10) 

where (p denotes the Fourier transform of ip. Quasi-neutrality, equation |2.6[ in Fourier 
space becomes 



^(k) = m J vdvJo{kv)g{k,v), (2.11) 

where 

= T, r-fTTTT- (2-12) 

(i + T)-ro(fc) 

Equations |2.5| and |2.6| are the starting point for all further investigations of this paper. 
This two-dimensional gyrokinetic system is the simplest paradigm for kinetic turbulence 
of magnetised plasmas. 



3. Collisionless invariants 



We focus on two quadratic invariants of the two-dimensional gyrokinetic system. To 
derive the first, let us start by noting that the average of g"^ is collisionlessly conserved 
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for each individual value of v. That is, the volume-averaged g^, as a function of velocity, 
must remain unchanged in the absence of collisions or forcing. Dividing by an arbitrary 
function of velocity, k{v), and integrating over v, a class of invariants is formed, which we 
will refer to as the "generalized free energy" (after this section, we use the abbreviation 
"free energy"): 

„y»,„)l = = 2./.<i„/^£||, ,3.1) 

where V is the system volume. A scaling theory for two-dimensional gyrokinetic turbu- 
lence may be developed, while allowing for arbitrary k{v) — under the constraint that 
k{v) varies only on the thermal velocity scale uth = 1 and so does not introduce ad- 
ditional velocity-space scales into the problem. Two functions, k = 1 and k = Fq, are 
particularly useful, so we introduce the following notation: 



Wgi=Wg[l], (3.2a) 
Wgo = Wg[Fo]. {3.2b) 



The second invariant of the gyrokinetic system, equations |2.5| and |2.6[ is particular to 
the electrostatic two-dimensional case: 



E : 



[{1 + T)if^ - ipToif] 



(3.3) 



1 r d^r 

2 J ~V 

We refer to this quantity as the "electrostatic invariant." Note that the operator Fq 
renders the second term negligible at large k (because in Fourier space Fq ^ k^^i see 
equation 2.9 1. For this and other reasons, it is convenient to define the quantity E(p 



(which is not in general an invariant): 



V 



(3.4) 



From equation 3.3 E{k) « (H- T)Etp{k) for fc > 1. 



We now demonstrate the conservation of Wg and E by deriving the "global budget" 



equations from the gyrokinetic system. Multiplying the gyrokinetic equation 2.5 by g/ n, 
averaging over gyro-centre position R and integrating over velocity, we get 



dt 



(3.5) 



where the injection rate of free energy is 

d^r 



I A ^^^^ 
It: / vdv 

V .] K 



(3.6) 

Multiplying the gyrokinetic equation |2.5| by {(p)^ and integrating over real-space co- 
ordinate r and velocity v, we obtain the global budget equation for the electrostatic 
invariant: 



where the injection rate of the electrostatic invariant is 



(3.7) 
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Sk = 



If 271 vdv 



(3.8) 



1 r d\ 

2J 'V 

It is useful to compare these results with what is known for three-dimensional gyroki- 
netics. Note that the electrostatic invariant is not conserved by the three-dimensional 

because of the particle streaming term v\\dh/dz that must be 



2.2 



gyrokinetic equation 

retained (just as the three-dimensional Navier-Stokes equation does not conserve enstro- 
phy, while in two dimensions it does) . Three-dimensional gyrokinetics does have an invari- 
ant that plays the role of the Navier-Stokes energy. The invariant is the perturbed free 
energy, W (see [Krommes fc Hu|[T994| Sugama et aL||1996 Howes et aL||2006 Schekochi- 



hin et al. 20081 2009). Its conservation is demonstrated by multiplying the gyrokinetic 



equation |2.2| by /i/Fq and integrating over velocity and position r. In our notation, W is 



W^Wgo + E. (3.9) 

For scales much smaller than the Larmor radius, we will see that the electrostatic invari- 
ant contributes negligibly to W and it is therefore justified to refer to either W or WgQ 
as the free energy. 



4. Charney Hasegawa Mima and Euler equations 

The CHM and Euler equations are popular paradigms for two-dimensional fluid tur- 
bulence. They can be derived as rigourous (and physical) asymptotic limits of the gy- 
rokinetic equation, as we now demonstrate. 

We begin with the gyrokinetic system [23| and [2?6l and take the long- wavelength limit, 
fc ^ 1, and the small temperature ratio limit r <C ij^ Formally, we assume fc^ ^ r and 
dt ~ C[.] ~ k'^Lp. The gyro-average in equation 
order fc^ in equation 2.6 



2.5 



acts as unity, while wc must keep Fq to 
(the zeroth-order part cancels), so Fq ~ V'^. The gyrokinetic 



system in this limit is 



dt 



vbo ■ V.g = C[g\ + J", 



27r / vdvg — Tip — V 



(4.1a) 
(4.1&) 

where V£;o = z x V(/j is the E x B velocity (without the gyro-average) . Now, by integrating 
equation [4.1 a over velocity space (noting that the integral of the collision operator is zero 
by particle conservation) and using quasi-neutrality, equation |4.1& we obtain a single 
equation for the electrostatic potential: 



dt(T - V^)^ + • V(-VV) = -FcHM, (4.2) 
where T^^^i = 2tt J vdvJ-. This is the inviscid CHM equation. The scale given by r^^/^ 
corresponds to the Rossby deformation radius in quasi-geostrophic turbulence or the ion 
sound radius, ps, in a plasma. The two-dimensional Euler equation is obtained in the 
case T = (the no- response model; see appendix [b| equation B 3|. 

There arc two invariants of the CHM/Euler equation that arc relevant to our discussion. 



The limit of small temperature ratio is somewhat restrictive but the CHM equation follows 
from it in a simple and transparent manner. We note that a more detailed system of fluid 
equations can be worked out for r of order unity but, for our purposes, it is not substantially 
diff'erent from the familiar CHM equation. 
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These are referred to as energy and enstrophy, although their physical interpretation 
depends on the specific scale of interest. They are found by multiplying equation |4.2| by 
(f and V^(/9, respectively, and integrating over the system volume: 



i^cHM "^(V + IVy'l'), (4.3a) 

Zc.M = \j ^(r|V<^p + (W)'). (4.3&) 
We now show that there is an additional invariant associated with the gyrokinetic de- 



scription (given by equation 4.1). Without loss of generality, we introduce the following 
ansatz for g: 



g^FoiT~V^)ip + g, 



(4.4) 



Quasi-neutrality (equation 4.16) implies that g does not contribute to particle density: 



vdv g = 0. 



(4.5) 



Now, multiplying the CHM equation 4.2 by Fq and subtracting it from equation |4.1a| 
we obtain a separate equation for g: 



dg 
dt 



- VBo • Vg = {C[g])^ + (J" - FoTcni. 



(4.6) 



Thus, g behaves as a passive scalar in the CHM/Euler limit, advected by the E x B 
flow determined by equation |4.2[ From equation |4.6[ it follows that there is another 
collisionless invariant: 



vdv 



k{v)' 



(4.7) 



where k{v) is an arbitrary function of velocity. The free energy invariant Wgo (equation 
3.26) can now be written in terms of the three gyrokinetic-CHM invariants: 



(4.8) 



'go = "•'ffO 

where we define WgQ = Wg[i^o]- As we have shown, each term here is conserved individ- 
ually in the CHM/Euler hmit. 

To complete our account of how the gyrokinetic invariants reduce to the CHM/Euler 
invariants in the long-wavelength limit, we note that the (gyrokinetic) electrostatic in- 
variant E reduces to the energy invariant i?cHM- From equation 3.3 taking Fq ~ 1 + V^, 
we have 



E = Ea^M. (4.9) 

As a final comment in this section, note that viscosity did not appear in the fiuid 
limit derived above. To recover the viscous CHM and the two-dimensional Navier-Stokes 
equation, one must raise the ordering of the coUisional operator. Specifically, the viscous 
term can be obtained by ordering C[.] ~ k^^dt- 
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5. Phenomenology 

In sections [6] [7] and |8j we will formally derive properties of fully developed gyrokinetic 
turbulence. In this section, we give a simpler derivation of some of the scaling results via 
a set of phenomenological arguments in the style of Obukhov ( 1941a|& ). We first sketch 
the phenomenology for stationary, driven CHM turbulence, and then outline how these 
arguments can be adapted for the essentially kinetic short-wavelength regime. 



5.1. GHM range, A; < 1 

Inspection of equation |4.3| reveals that the two invariants of the CHM and Euler equations 
have spectra that are constrained to satisfy Zcmmik) / Ecmmik) = k^ ■ By the argument 



due to Fj0rtoft (1953), this implies that a dual cascade will occur. That is, injection at 
a scale ki will give rise to two cascades: an inverse cascade of energy to k < ki and a 
forward cascade of enstrophy to k > ki. 

Using the assumptions of isotropy and local (in scale) interactions, the CHM equation 
4.2 is translated into phenomenological language as follows: 



u^^T + r^)iff, ^ ip^,r\ (5.1) 

where w^l is the nonlinear decorrelation rate and ipi is the characteristic amplitude of 
(fi at scale £ ~ k~^. The essential assumption that leads to a scaling theory is that the 
flux of enstrophy is constant for k > fcj , while the flux of energy is constant for k < ki 
Multiplying the right hand side of equation 



5.1 



by i'^ (pg gives the nonlinear flux of 
enstrophy, which must be equal to the total rate of enstrophy injection. This flux must 
be independent of scale in the stationary state: ^ £z = const. Therefore, 



,1/3.2 



(5.2) 

and the spectrum of potential fluctuations]^ is E^{k) ^ fc~^. 

To flnd the scaling law for the inverse-cascade range, k < ki, we multiply the right- 
by ipf to get the energy flux, and set this equal to the energy 



hand side of equation 
injection: i 



5.1 



const. This implies 



-l/3«4/3 



(5.3) 



and, therefore, Em(k) ^ k ^^/^ 



Finally, the cascade of Wg (equation 4.7 1 is passive and directed forward. Constancy 



of nonlinear flux of this invariant implies t 
have for k > h 



-2-2 

91 'Pe 



e — const Using equation 



5.2 



we 



9i 



-1/2 - 

e ' e 



-1/6 



const. 



(5.4) 



Therefore, the spectrum of fluctuations of g scales as Wg(k) ~ A: , the usual Batch- 



elor (1959) scaling for a passive scalar (also, see Lesieur & Herring 1985 for the two- 



dimensional case). 



^ The spectra -EcHM(fc) and Zci\m(}i) have different scaling laws depending on the relative 
sizes of I ? an d r. Thus, it is more convenient to speak about the spectrum E(p{k) (defined by 
equation 3.4 1, which has a single scaling law. 

* Note tEat the passive cascade at A; > fc^ i s m arginally non-local because the decorrelation 



rate scales as cot^L ~ 'fiel^ , which, by equation 5.2 implies that all scales in the forward-cascade 



range contribute similarly to the time evolution of g. 
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Figure 2. Correlation in velocity space. Contours of constant potential are plotted in a blue 
color scheme. Overlaid are the orbits of two particles that have the same gyro-centre. During a 
nonlinear turnover time, the phase- mixing process will de-correlate the motion of these particles 
because their Larmor radii differ by the correlation length £ of the potential. Thus, ~ £. 
(Colour online.) 



5.2. Nonlinear phase-mixing range, fc ^ 1 



Schekochihin et al. (2008 2009) proposed a phenomenological description of the large- 



k forward cascade in gyrokinetic turbulence. We reproduce their arguments here, in a 
slightly different form, and also propose a phenomenology to describe the inverse cascade. 
This will motivate the formal derivations of sections [6] and [8] and highlight the saUent 
physical processes. 

In gyrokinetics, the drift motion of particles is determined by fluctuating fields aver- 
aged along the particle Larmor orbits. Two particles sharing the same gyro-centre but 
having different velocities perpendicular to the equilibrium magnetic field are subject to 
a different averaged electrostatic potential and thus have distinct E x B drifts. This is 
illustrated in figure [2j One can see from this cartoon that if the E x B flow has a cor- 
relation length £, then the particle motion must be decorrelated over the velocity-space 
scale £v such that £v ^ I. 

Now let gi(v) be the characteristic amplitude of the distribution function g at scale 
£. In view of the above argument, we assume this quantity to have random structure in 
velocity space on the scale ~ £. From quasi-neutrality (equation 2.6 ) we can estimate 
the scaling of tp^ in terms of gi. We make the substitution k^^ £ in the Fourier-space 
expression 2.11 and take fc ^ 1 (so that /3(fc) is a constant). The relationship between 



and gi is then derived as follows: 
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vdvJo{v/e)ge{v) ^ ^1/2 / v^/^dv cos{v/e - 7r/4)ge{v) ^ e^^^el^^ge, (5.5) 



where we have used the large-argument approximation of the Bessel function Jq. Be- 
cause of the random velocity dependence of ge{v), we argue that the velocity integral 
accumulates like a random walk — the number of steps scales as and the step-size 

1 /2 

scales clS £y , SO the typical "displacement" scales as £v ■ Note that in the final expres- 
sion we have introduced the amplitude gi (distinct from ge{v)), which can roughly be 
interpreted as the root- mean-square value of the fluctuations in gi(v) over velocity space: 

gi ^ Jvdvgj{v). Taking iy ^ i (as motivated by the arguments illustrated in figure |2| 
in equation 1 5. 5 [ we find 



^gi- (5.6) 

This immediately implies that there is a possibility of a dual cascade for the two 

h olds, the spectr a Wg{k) and E{k) must satisfy 



invariants Wg and E. Indeed, if 



5.6 



Wg(fc) ~ k'^E{k). By the argument of Fj0rtoft (1953), this constraint leads to a dual 



cascade where E is transferred to larger scales and Wg is transferred to smaller scales. 
Thus, if the injection scale is ki 1, we predict the existence of two cascades: a forward 
cascade of Wg to k > ki and an inverse cascade of to A; < fc,;. 



From the gyrokinetic equation |2.5( we estimate the nonlinear decorrelation rate from 
the nonlinear term by substituting • V — t- Jq{v/£)(Pi£~'^ , where the Bessel function 



comes from the gyro-average in equation 2.10 Then the nonlinear decorrelation rate is 



c^«L - r^^^'Pe, (5.7) 
where we have used the large argument approximation for Jo{v/£) to estimate that the 
gyro-averaging introduces a factor of i^^^. To determine scalings in the forward-cascade 
range, we assume locality of interactions and constancy of the flux of free energy (see 
section [6j equation 6.11 for the exact statement of constancy of flux). This gives 



2 



£ ^^^^{gj ~ = const. 
Using equation 5.6 we have gft~^/'^ ~ and, therefore. 



(5.8) 



"Pi 



^l/3„l/6 

^ sT£'l'. 



(2009 2010&) 



8.2 



(5.9) 
(5.10) 
these are 



The resulting spectra are Wg{k) cx k ^/'^ and E{k) cx k ^^/-^ (in section 
derived more formally). These spectra have been confirmed numerically by Tatsuno et al. 



The analysis of the inverse-cascade range, k < ki, follows as in section [5TT| Assuming 
constancy of the fiux of the electrostatic invariant E in the inverse-cascade range (see 
section 



6.4 equation 6.18), we have, using equation 5.7 



Hence 



£e- 



(5.11) 



.l/3«l/2 



(5.12) 
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The corresponding spectrum is E{k) ~ Etp{k) ^ k ^. If we assume that the relationship 
5.6| also holds in the inverse-cascade range, equation |5.12| implies 



9i 



,l/3«-l/2 



and so the free-energy spectrum in this range is Wg{k) 
formal derivation of these scalings. 



See section 



8.4 



(5.13) 
for a more 



We note that the numerical experiments of Tatsuno et al. (2010a I have confirmed that 
the locality of the cascade is a good assumption and that the relationship ~ £ is well 
satisfied for both the forward and inverse cascades. 

5.3. Dissipation scale 

Having derived the inertial-range scalings for the forward cascade, it is easy to estimate 
the collisional cutoff scale that marks the end of the inertial range. For this purpose, 
we notice that the collision operator acts by second derivatives in phase space, and so 
assuming £y ~ i, we may take the collisional damping rate at scale £ to be Wcoii ^ v£~'^ , 
where v is the collision frequency. At the collisional cutoff scale tc, the nonlinear and 
dissipation rates must balance: Wnl(^c 
obtain 



ajcoii(^c)- Using equations 5.7 and 5.10 



we 



1/5 



(5.14) 



which is the gyrokinetic analog of the Kolmogorov scale. Note that the cutoff scale may 
be written in terms of Do, the Borland number ( Tatsuno et aL||2009 ). Setting i — \ va. 



equation 



5.10 



one finds ~ v\ ^^^d so kc ~ Do^^^, where Do = ^pi/v- Do plays the role 
of the Reynolds number for the nonlinear phase-space cascade in gyrokinetic turbulence, 
where the outer scale of the turbulence is taken to be the Larmor scale because it marks 
the start of the nonlinear phase-mixing rangej^ 

5.4. Gyrokinetic cascade through multiple scales 

Using the phenomenological scalings derived above, we now give a sketch of the fully- 
developed cascade from the low-fc CHM/Euler limit to the high-fc nonlinear phase-mixing 
limit. Note that the kinetic invariants Wg and E are exact invariants for both the high-fc 
and low-fc regimes, so they cascade through this entire range, including fc ^ 1 (there is 
no dissipation at the Larmor scale). We may use what we have learned about the rela- 
tionship between the CHM/Euler equations and our two-dimensional gyrokinetic system 
to describe this cascade in detail. In particular, equation |4.8| states that the gyrokinetic 
invariant Wg is composed of three separately conserved quantities in the CHM limit. In 
spectral form, the equation is 

Wg{k) = Wg{k) + TEc^^,ik) + ZcHM(fc) = Wg{k) + [t^ + 2t k'' + k^] E^{k). (5.15) 
For the electrostatic invariant, equation |4.9| gives 



E{k) « Ea^Uk) - (r + k^)E^{k). (5.16) 
The cascades traverse two special physical scales, the ion Larmor radius (fc — 1) and 



' In dimensional units, Do — ujkl{p)/l' = ipp/ Bop^v, where ajNL(p) is the nonlinear decorre- 
lation rate and ipp is the fluctuation amplitude, both evaluated at the Larmor scale. 
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the ion sound radius {k — r^/^). While nothing is dissipated at these scales, the physics 
of the turbulent fluctuations changes, giving rise to different scaling regimes. 

The range k <C t^/^ is referred to as the 'potential limit' of the CHM equation, because 
the spectrum of i?cHM becomes proportional to the potential fluctuation spectrum, E^[k). 
The range r^/^ ^ fc ^ 1 is called the 'kinetic limit' of the CHM equation because Ecmm 
may be interpreted as the kinetic energy density of the flow. The range 1 <^ k <^ kc 



(where k^. is the dissipation scale; see section 5.3 ) defines what we have been referring to 
as the nonlinear phase-mixing regime. 

We have assumed that there is localised injection of free energy at a single wavenumber 
ki. We now examine two possible scenarios: ki <C r^/^ (injection at large scales) and 
fci 1 (injection at small scales), depicted in figures [3] and [4j respectively. 



5.4.1. Injection at large scales 

In figure [3| the injection is at large scales, in the potential limit of the CHM equation, 
where ki <C r^^^. In this limit, the electrostatic invariant E reduces to E k tEi^. An 

with the scaling derived 



inverse cascade of the electrostatic invariant occurs at fc <C A:^ 

-11/3 



5.1 



(equation 5.3): E{k) - Em{k) - k 



in the section 

For fcj ^ fc ^ 1, a forward cascade of the free energy occurs, which takes the form 
of the cascade of enstrophy Z and the passive cascade of Wg (the invariant due to 

The spectrum of the potential fluctuations 

is Eu 



(fc) 



the zero-density part of g, equ ation 4.4). 

(equation 5.2). To determine the scahng of Wg(fc), we must assume 



something about the strength at which Wg is driven. Let us suppose that the forcing 



drives Wg weakly enough, such that we may neglect Wg{k) in equation 5.15 
for k < ri/^ we have Wg{k) « TE{k) « T^E^{k) and for r^/^ 



Then, 

^ A: ^ 1 we have 



Wg{k) « Zcmm{k) ~ k'^E^{k). The forward cascades of Zchm and Wg will combine 
around Larmor radius scale k ^ 1, to drive the cascade of Wg in the nonlinear phase- 
mixing regime k ^ 1, giving Wg{k) ~ fc"'*/'^ and E(k) ^ Eip{k) ^ k~-^^^^ (see section 



5.2 equations 5.9 and 5.10). These scahngs are summarised in figure |3] 



5.4.2. Injection at small scales 

In figure |4j the injection is at small scales, in the nonlinear phase-mixing range, fc 3> 1 



The free energy Wg cascades forward to the coUisional cutoff scale, kc (equation 5.14), 



whereas the electrostatic invariant E cascades inversely. For 1 ^ A; <C fcj, the scaling of 



E{k ) and Wg{k) are taken directly from the phenomenology (section 5.2 equations 5.12 
Eik) ~ k-^ and Wg{k) 



and 



5.13) 



k'> 



Again, the identity of Wg changes in the long wavelength limit, fc <C 1. All fluctuations 
at these scales are due to the inverse cascade of E, so Wg{k) is determined from Eip(k) 
The inverse cascade is characterised by the scaling E(n{k) ~ 



by via equation 
(equation |5.3 



5.15 



and E{k) « k^Eipik) - fc-^/^. In 
T^E^ik) - fc^"/3^ E{k) « TEu,ik) 



We argue that there is no transfer of Wg from small scales because g is 
passive. Thus Wg{k) is absent in equation 5.15 In the kinetic limit of the CHM equation, 
r < A:^ < 1, we have Wg sa Zchm(A:) ~ k*E^{k) 
the potential limit fc^ <C r, we have Wg w TE{k) 
k ' and .Z^chm 

(A:) « Tk'^Eip{k) k^^/^. These scahngs are summarised in figureji] 
As we have shown in this section, much about the gyrokinetic turbulence can be 
understood by simple heuristic arguments. Yet it is important to put these results onto 
firm foundation by making more precise, quantitative derivations. We also would like to 
develop a more detailed theory of the phase-space structure of turbulent fluctuations. 
This is done in sections |6] [7] and [8j 




6. Statistical theory 

In this section, we will analyse the gyrokinetic system, equations 2.5 and |2.6[ using 
familiar statistical methods from fluid turbulence theory. In sections 6.1 and |6.2[ we dis- 
cuss symmetries of gyrokinetic system and the related statistical assumptions. Sections 
|6.3| and |6.4| will be devoted to deriving exact statistical scaling laws in analogy to Kol- 
mogorov's 'four-fifths' law (or the 'three-halves' law for two-dimensional NS turbulence). 
There will be two scaling laws: one for the forward cascade of the free energy Wg and 
one for the inverse cascade of the electrostatic invariant E. Note that the exact laws 
constitute a precise statement of the constancy of nonlinear flux, which was used in the 
phenomenology of section |5.2| 



6.1. Symmetries of the gyrokinetic system 
Statistical arguments in turbulence theory are anchored to the underlying symmetries 



of the dynamical equations of interest (see Frisch 1995). Although specific boundary 
conditions, initial conditions and forcing mechanisms break these symmetries at large 
scales, a guiding principle of turbulence theory is that at sufficiently small scales, far 
from boundaries, a fully developed state will emerge where the symmetries of the dy- 
namical equations are restored in a statistical sense. We will appeal to this principle in 
the arguments of this paper. 

As in fluid turbulence, one uses the translational invariance in time and position 
space for arguing statistical homogeneity and stationarity respectively; likewise, statisti- 
cal isotropy originates from the rotation and reflection symmetries. For gyrokinetics, the 
scaling symmetry is also particularly interesting. In the coUisionless limit, and for k ^ 1, 
the gyrokinetic system is invariant under the transformation 
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Figure 4. Scenario 2: Injection at small scales. The scaling indices for each sub-range are 
printed in bubbles over the plots of the invariants. (Colour online.) 



[g, V, I", v) (g, /LfV, fJ'V), (6.1) 
where /i is a scaling factor. We will ar gue t hat this symmetry implies a dual-scaling in 
velocity and position space (see section 8.1 ). p°] 



6.2. The ensemble average 

We denote the ensemble average with an over-bar and define the correlation function of 
g between two points in phase space, (R, v) and (R', w'), as 



G = gill,v)gill',v')=gg'. (6.2) 
It is also useful to define the second order structure function: 

^2 = V, (6.3) 
where Sg = g(R', v')^g{'R, v). Homogeneity in position space (gyro-centre space) implies 
that G and S2 are functions only of the increment = R' — R. In contrast, velocity 
space is fundamentally inhomogeneous, as there is no translational invariance in itj^ 
However, we expect approximate "local homogeneity" in velocity space for sufficiently 



This symmetry also provides a formal basis for the phenomenological argument £v ~ £. 
See section 5.2 figure [2] 

Velocity dependence in the Maxwellian and the gyro-average both break translational 
invariance. 
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small velocity increments £y — v' — v. Specifically, the two-point statistical quantities, 
such as structure and correlation functions, should vary strongly with the increment iy 
but weakly (smoothly) with the centred velocity v — {v' + v)/2. 

6.3. Third-order Kolmogorov relation for free energy 
We first derive the budget equation at scales £ and We take the gyrokinetic equation 



2.5 for g and g' and multiply by g' and g, respectively, then add the two resulting 



equations to obtain 

f)r< 1 

— - - V, • = g {Cmw + .9' (CW>K + .9-^' + .9'-^, (6.4) 

where the third-order structure function is defined 



S'<"' = 6vESg^ (6.5) 



and — — v^;. To get equation 6.4 we have used incompressibility of v^; (V • v^; 



0) and statistical homogeneity to manipulate the nonlinear flux to obtain — g'v e ■ Vg + 

9^'e ■ V'5' = V, • 5vbV/2. 

Now if we take £ = 0, — 0, the nonlinear term in equation |6.4| is zero by homogeneity. 
The result is the kinetic 'global budget' equation: 

dg 



^2g{C[h])^ + 2gT. (6.6) 

Assuming stationarity, the left-hand side of the equation is zero and we have the steady 
state balance between forcing and dissipation 



gT= -g{C[h])^ = e,Av), (6.7) 
where we have defined the kinetic free energy injection rate e„/(w). 

We now assume that the forcing is restricted to a narrow range around an injection 
scale which we denote £i. For £ <^ ii, <^ £i and ^t, <C 1, the forcing term in the budget 



equation 6.4 reduces to 



gJ^' + g'J'^2gJ^^2e„{v). (6.8) 

Now we take the limit of small collisionality: the collision operator is proportional to v, 
the collision frequency, which we take to b e sm all so that, for fixed £ and £y, the coUisional 
term is negligible in the budget equation 6.4 ^ Thus, by assuming small collisionality, 



statistical stationarity, £y ^ 1 and £ £i we have from equation 6.4 



Ve-S':^'' ^-Ae„{v). (6.9) 
Note that although homogeneity and stationarity have been assumed, we have not yet 
assumed isotropy. Assuming it now, we may write S^^-' = tS'"^'\£), where 1 = £/£. Then 



equation 6.9 reduces to scalar form: 

1 d 



^^^{£Sr)--^eAv). (6.10) 



Wliat is implicit to this assumption is that a collisional cutofT scale £c exists, such that 
for £,lv ^ Ic, the collisional term may be neglected — the existence of this scale can either be 
demonst rated a posteriori or assumed from the arguments of the phenomenology section (see 
equation 5.141. 
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The solution to equation 



6.10 



5<»-) = -2e„{v)L 



(6.11) 



6.4. Third-order Kolmogorov relation for the electrostatic invariant 

To derive the scale-by-scale budget equation for the electrostatic invariant E we proceed 
similarly to the analysis of section [673l However, in this case, we are dealing with functions 
of real-space variable r and its increment tr = v' — r. Also, velocity will appear only as 
a variable of integration, so it will simplify our notation to refer to a single variable v in 
what follows — i.e. R = r — p{v), R' = r' — p{v) and g' = g(R', v), etc., from which it 
follows that i = and we will refer to a single increment in position space, t. 



We take the gyrokinetic equation 2.5 for g and g' and multiply by ip' and ip respectively, 
then ensemble-average and integrate over velocity v. The result is 



(i + r-ro)^-v,.5r' = 



27r ^' / vdv + 27r / vdv + + p'T^, (6.12) 



where we have defined the correlation function $ = pip' and introduced the third-order 
structure function 



(f) — 2n vdv 5 {Lp)-^SvESg, 



(6.13) 



where 6 (</5)r = {'p')r' — {'p)r- Note that by homogeneity, we have Fo^" = T'q^ 
where Fg operates in ^-space. Lastly, we have defined the electrostatic forcing Te 
277 J vdv (T)^. 



We can express the first term on the right side of equation 6.12 in terms of the rate 
of change of the invariant E and structure functions involving ip. Using homogeneity, we 
find the following relation 



(l + r-Fo)$ = 2S 



Sip5{TQip) 



(6.14) 



where S{Tqip) = T^ip' — Toip. In the theory of fluid turbulence, it is standard at this 
stage to neglect the time derivatives of structure functions in the stationary limit {e.g. 
Kolmogorov||194f ) , assuming essentially that a self-similar regime is attained. We follow 
this practice and substitute {1 + t — Fo)9t<f> w 2dtE into equation 6.12 Taking the limit 



of small collisionality, equation |6.12 becomes 



dE 
~dt 



1. 



■^3 



1 



1- 



(6.15) 



2 ^ 2" 2' 

For £ — 0, we have the global balance of E, which describes the continual accumulation 
due to forcing: 



dE 
~dt 



ipJ^E 



(6.16) 



where Ee is the injection rate of the electrostatic invariant E. Recall that the forcing 
is assumed to be restricted to a scale £i as explained in the section |6.3[ For the inverse 
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cascade, we are considering an incrtial range such that £ ^ £i. In this Umit, the forcing 



term in equation 6.15 is zero and we find the following exact (anisotropic) result: 



V, • S 



- 2e. 



(6.17) 



Assuming isotropy, we have S^^^ — £S^'^\i) and so equation 6.17 gives 



^3 



(6.18) 

Note that in the CHM/Euler limit (£ ^ 1), equation 6.17 may be manipulated to coincide 



with t he result of Boffetta et al. { 2002 ) for CHM turbulence and the result of Bernard 
(1999) for two-dimensional NS turbulence: •(5u|(5uP — 4e^, where u = veo = zx 'Vip. 



The two exact laws derived in this section, equations |6.11| and equation |6.18[ hold 
for arbitrary spatial scales, i.e., they hold above, below and at the Larmor radius scale. 



(Note, however, that for equation 6.11 we have assumed that the velocity scale satisfies 
£v <C 1.) In the sections that follow, we will consider the limit £,£v <C 1. This is the 



essentially kinetic limit (see section 5.2 for phenomenology) where nonlinear phase-mixing 



is dominant and the turbulent cascade involves both velocity and position space. 



7. Spectral formalism 

In this section we introduce a way of characterising velocity scales using a zeroth-order 
Hankel transform. This produces an elegant spectral representation of gyrokinetics, where 
the electrostatic potential and electrostatic invariant are compactly expressed in terms of 
the Hankel transform of the distribution function. In section [8] we will use this formalism 
to develop a spectral scaling theory for the inertial ranges in phase-space. 

7.1. Hankel transform 
The Hankel transform of a function g{v) is 

/•oo 

dip) = / vdvJo{pv)g{v). (7.1) 







The Hankel transform is its own inverse. This is easily proved using the orthogonality of 
the zeroth-order Bessel functions: 



Ja{pix)Jo{p2x)xdx = S{pi -p2)/pi. (7.2) 
Using this identity, the Plancherel theorem for two functions / and g is 







vdvf{v)g{v) = / pdpf{p)g{p). (7.3) 
Jo 

The following integral involving three Bessel functions will be useful in deriving the 
mode coupling due to the nonlinearity: 



oo 



vdvjQ{piv)Joip2v)JQip3v) ^ K{pi,p2,P3), (7.4) 







where if pi, p2 and p^ do not form the sides of a triangle, then = 0; if they do, then 
K = l/(27rA), where A is the area of that triangle. 
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7.2. Spectral gyrokinetic equation 

With the help of the above identities, we may now rewrite the gyrokinetic system as a sin- 
gle equation in Hankel-Fourier space. The Hankel-Fourier transform of the distribution 
function (/(R, w) is 

g{k,p) = ^ f d^n f vdvJo{pv)e-'''-^g{K,v). (7.5) 
Jr Jo 

Quasi-neutrality, equation |2.11[ can now be expressed simply as 

^(k)=/3(k)g(k,fc). (7.6) 
Neglecting collisions and applying the Hankel-Fourier transform to the gyrokinetic equa- 



tion 2.5 we find the spectral form of the gyrokinetic system 
dg{k,p) 



dt 2tt 



^ J d^k'Pik') k • (z X k') y qdq K{k\p, q) g{k! , k')g{k - k', q), (7.7) 



where fi is defined in equation 2.12 and K is given in equation 7.4 



7.3. Free energy spectrum 

Let us apply this formalism to define the spectral density of free energy. First, we use 
homogeneity and isotropy to express the conserved quantity g'^ in terms of the correlation 
function G{t, v' , v): 

'ffiv) = J kdkpdpidh' dv' Ja{ki)Jo{pv')Jo{pv)G{(i, v' , u), (7.8) 



where we have used equation 7.2 and the identity Jd'^ke*'" = (27r)^(5(r). Note that 
isotropy has been used to reduce the Fourier transform in vector position increment ^ 
to a Hankel transform in scalar increment £. For simplicity, we consider generalised free 



energy corresponding to k = 1 (see equation 3.2a 



/2 
vdv^-. (7.9) 

However, we stress that the spectral scaling results derived below will hold for general 
k(v) as long as this function has weak dependence on velocity {i.e. for k{v) with ^„ ~ 1). 
Thus, scahng results will be reported in terms of Wg{k,p), although Wgi will be used 
during intermediate steps, for the sake of consistent notation. 



From equation 7.8 it follows that 

Wgl = j dkdpWgl{k,p), (7.10) 



where 



Wgi{k,p) = TTpk J £devdvv'dv'Jo{ke)Jo{pv')Jo{pv)G{i,v',v). (7.11) 
Wgi{k,p) can be re-expressed in more convenient velocity variables: the increment vari- 



able iy = v' — V and the centred velocity v = {v + v')/2. Then equation 7.11 becomes 



Wgi{k,p) = wpk / £di{v^ - fjA)deydvJo{k£)Q{p,£y,v)G{£,£y,v), (7.12) 
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where Q{p,L,v) = Jo[p(2w + £„)/2] Jo[p(2w - 4)/2]. 

Note that equation |7.11| is not invertible, unHke the case of the Wiencr-Khinchin for- 
mula in fluid turbulence — i.e., the correlation function G{i,v',v) cannot be recovered 
from Wg{k,p). Because of the lack of homogeneity in v, we have integrated over v and 
v' to define scales in velocity space via a single velocity-space "wavenumber" p. 

The spectral density can also be expressed in terms of the Hankel-Fourier transform 



of the distribution function (using equations 7.3 and 7.5) 



Wgiik,p) ^TTpk g{k,pY. 
7.4. Spectrum of the electrostatic invariant 



(7.13) 



Using equations 3.3 and 7.6 the spectral density of the electrostatic invariant can now 
be expressed 



Eik) - ^^mW ^ ^Wg^{k,k), 



(7.14) 



which also implies that E^{k) = /3'^{k)Wgi{k,k)/{2Trk). We can also write S2E = I^veI"^ 
in terms of the spectral free energy density Wgi{k,p). Using equation 



7.6 



we find 



S2E^ J dkH{k,i,v,v')Wgi{k,k), (7.15) 
where H{k,e,v,v') = -kl3^{k)/TT [JUkv') + J§{kv)~2Ja{kv')Jo{kv)Jo{ke)]. 



8. Spectral scaling theory 

In sectionjTj we did not make any assumption about the size of fc or p and the equations 
presented were exact. In this section, we will derive some approximate spectral scaling 
results at the sub-Larmor scales, k ^ 1 and p ^ 1. In section |8.1| we introduce a 
similarity hypothesis for phase space that will allow us to infer spectral scalings from the 
exact laws of section [6] In section [8?2{ we derive general relationships between the scaling 
of the structure function, 5*2, and spectral densities Wg{k) and E{k). Then, in section 
|8.3[ we derive specific scaling laws for the forward-cascade range and, in section [8^ for 
the inverse-cascade range. 

8.1. Self- similarity hypothesis and scaling in velocity space 

In the theory of fluid turbulence, Kolmolgorov's four-fifths law can be used to infer the 
scaling index for the velocity field under the assumption of self-similarity. This allows 
derivation of the (approximate) energy spectrum directly from the exact result for the 
third-order structure function. We will follow an analogous approach for gyrokinetics in 
the following sections. In our case, however, we would like to determine scaling with 
respect to the increment in addition to i (and corresponding phase-space spectral 



scaling laws). To accomplish this, we supplement the exact laws, equations 6.11 and 



6.18[ with the following dual self- similarity hypothesis: given fixed scales £ and £v, such 



that £v,£ <^ 1, the increment Sg obeys 

Sg{X£,X£,)^X''^Sg{£,£,). (8.1) 



The notation = is meant to echo the "equality in law" of Frisch ( 1995 1. This means that 



6g exhibits this scaling statistically, i.e. whenever it appears inside an ensemble average. 
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8.1 



is true 



For instance, it implies S2{\£, \£v) / S2{£, iv) — A We argue that equation 
on symmetry grounds. In particular, since the coUisionless gyrokinetic system is invariant 



to simultaneous scaling of R and v by the same factor (this is the symmetry 6.1), the 
statistics of the stationary state should obey this symmetryf^ 
As we prove below (see equations 8.7 and 



the dual self-similarity hypothesis for 



dg implies self-similarity for the gyro-averaged E x B velocity difference Sve- 



.2. Spectral scaling in the inertial ranges 



(8.2) 



Consider equation |7.12| To begin, we take 1 so that we can use the large-argument 

approximation of the Bessel functions in Q{p, £v, v). Also, we assume that the correlation 
function G{£, ly, v) is peaked at small iy, so the integral is dominated by ^ v. Since the 
integral is dominated by small i^, the bounds of the integration over £y are unimportant. 
In this limit, 



1 



[sin{2pv) + cos{ply)] 



(8.3) 



Substituting this expression into equation |7.12[ we see that the term proportional to 
sin(2pz;) oscillates rapidly in v — while other quantities in the integrand depend only 
weakly on v — so this term is neghgiblep^ Then, using G = [g'^ + {g'Y ~ 'S'2]/2, we can 
express Wg{k,p) in terms of 5*2 = Sg"^, for non-zero k or p: 



Wgi{k,p) w / ld£ vdv d£yJo{k£)cos{p£y) S2{£,£v,v) 



(8.4) 



Then, dual self-similarity of 5g (equation 8.1) implies that Wgi{k,p)^ and, therefore, 
Wg{k,p), must satisfy the scaling 



Wg{\k,\p) = X''Wg{k,p), 



.5) 



with 



v = -2-2hg. (8.6) 
We will now determine the relationship between the scaling indices v and (defined in 



equation 8.2 1. Expanding equation 7.15 in the limit few 3> 1 and £y <C w, as done for Q 



above, we get 



S2E 



dk [1 - cos{k£y)Jo{k£)] Wgi{k, k). 



•7) 



Using equations 8.7 and 8.5 we find 



Note that, although 5g must depend on the centred velocity its dependence is weak, and 
thus we argue that the scaling symmetry is only in I and 1^. The sensitive dependence in 1^ is 
due to the fact that the fluctuations must decorrelate in veloc ity s pace on the scale 1^ ~ t, as 



argued by Schekochihin et al.\ (20081 and reiterated in section 5.2 



with illustration in figure [2] 



One can think of this as "local homogeneity," though it is clearly difi'erent from homogeneity m 
the sense of neutral fluid turbulence. 

Thi s stat ement can be confirmed a posteriori by considering the scaling results given by 
equation 8.11 For instance, if we consider ^ = 0, we then have S2 ~ and the term we 

have neglected is indeed smaller by a factor of p~^. 
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= -1-v. (8.8) 

8.3. Spectral scaling in the forward-cascade range 
As a consequence of the dual self-similarity hypothesis and the exact third-order law, 



equation 6.11 we have the following scaling law for 5*3 



5<"''(A£, A4) = A5<"''(£,4) (8.9) 
Thus, the scaling indices for dg and Sve must satisfy the following relation in the forward- 
cascade range: 



2h„ 



Combining this with equations |8.6| and |8.8[ we have 



(8.10) 



And, therefore, taking p 



hg = 1/6, 
h^ = 2/3, 
v = -7/3. 

k, we have the power law 
Wg{k,k) cx fc-^/^ 



This implies (see equation 7.14) that in the large-fc limit 

1 



E{k) « (l + r)£;<^(fc) 



fc(l 



Wgi{k,k) cx fc'i°/3. 



(8.11a) 
(8.11&) 
(8.11c) 



.12) 



.13) 



This agrees with the E(p{k) scaling derived phenomenologically in section 5.2 (equation 



5.10 also, see Schekochihin et al. 


2008 


Tatsuno et al. ( 


2009 


20106). 



8.3.1. Phase-space spectrum Wg{k,p) 

We have found that the phase-space spectrum Wg{k,p) obeys the scaling law given 
by equation 8.5 with v = —7/3. To obtain this result, we have assumed k ^ 1 and 
p ^ 1, but have assumed nothing about the relative size of k and p. In the subsidiary 
limits k ^ p and k -^p, we can show that the spectrum Wg{k,p) has power law scalings 
separately in k and p. These scalings, obtained in appendix [Cl are 



Wg{k,p) cx 



for /c>]3>l, 
for p > fc > 1. 



(8.14) 



The one-dimensional spectra Wg{k) and Wg{p) {i.e., Wg{k,p) integrated over p and k, 
respectively) are also derived in appendix [C] We find the following power laws: 



Wg{k) CX 



.15) 



confirming the arguments of section 5.2 (equation 5.9), and 



Wg{p) cx p-^'\ 



(8.16) 
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The power law for Wg{k) agrees with the prediction by Schekochihin et at. (2008 2009), 
while the power law for W^(p) is a n ew prediction. Both scaling laws have been con- 



firmed numerically by Tatsuno et al. (2009 201061. As a consistency check one may 



integrate the asymptotic forms from equation 8.141 directly, to obtain roughly Wg{p) 



/jf dkk-^/^p-^ + J°° dkk-^p-^/^ ~ p-^/^ and similarly for Wg{k) 



8.4. Spectral scaling in ttie inverse-cascade range 
We now turn to the scaling of the spectra for the high-fc inverse-cascade range, 1 <C 



fc ^ ki. As in section 8.3 we will combine the dual self-similarity hypothesis with an 



exact law; in this case, the exact law is equation 6.18 the scaling law for Si'^\ However, 



there is now a complication arising from the velocity integration in the definition of S''^'' , 
equation |6.13| As we have argued in section |5.2[ integration of the distribution function 
effectively averages over the random velocity dependence and, thus, acts as a partial 
average, which we cannot assume commutes with the formal ensemble average. 

From the random- walk argument of section [5?2| we expect that the velocity integration 
of Sg should, roughly speaking, introduce a factor of i^^^. This 'rule of thumb' has been 



validated in the forward-cascade range by the numerical results of Tatsuno et al. ( 2009 



201061; furthermore, it agrees with the results obtained using the Hankel transform 



treatment. To formalise this, we introduce the scaling hypothesis 



dv F{v) dg{Xe, v) = A'^'+i/^ / dv F{v) 5g{l, v), 



.17) 



where F(v) is an arbitrary function, varying over scale ~ I. Also note that Sg{£,v) 
g{'R! ,v) — g{Il,v) depends on a single velocity v. 



Applying the hypothesis 8.17 to the result for 5*3 , equation 6.18 we find the following 



relationship between the scaling indices 



2/ie + hg = -1/2. 
Combing this with equations |8.6| and |8.8[ we get 



(8.18) 



-1. 



(8.19) 



Therefore 



Wg{k,k) cx k-\ 
E{k) cx k-^. 



(8.20a) 
(8.206) 



Equation |8 . 20 6| agrees with the phenomenological result of section [53] equation |5.12[ Let 
us now determine the free-energy spectral scalings. First note that by equations |8.6| and 
|8.19[ we have 



hg = -1/2. 



(8.21) 
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Now, if the free energy spectrum is a power law in k-p space as we have assumed in the 



forward cascade range (see equation 8.14), then we find 



Wg{k,p) cx 
Wg{k) cx /c°. 




k:^ p 

k<^p 



(8.22a) 

(8.226) 
(8.22c) 



where we have adapted the arguments of section 8.3.1 for the inverse-cascade result 
hn = —1/2, expanding equation 8.4 in the limits <C ^ and iy ^ £. Equation 8.226 



agrees with the phenomenological result of section |5.2[ equation |5.13[ 



9. Conclusion 

9.1. Summary 

In this paper, we have considered forced two-dimensional gyrokinetics as a simple paradigm 
for magnetised plasma turbulence. We have explored how the nonlinear phase-mixing 
mechanism gives rise to a phase-space cascade of free energy to scales much smaller than 
the Larmor radius. We have investigated the relationship between the inertial-range 
cascades in fluid theories {i.e. the CHM and two-dimensional Euler equations) and gy- 
rokinetics, and found a simple relationship between the fluid invariants and the kinetic 
invariants, given by equation 



4.8 



The gyrokinetic free energy invariant Wg (equation 
S.ll) is cascaded to fine scales and is ultimately dissipated by collisions. We predict that, 



concurrently, the gyrokinetic electrostatic invariant E (equation 3.3) cascades inversely 
from very fine scales {e.g., the electron Larmor radius) to long wavelengths — possibly 
even to wavelengths of the CHM/Eulcr limits, much larger than the ion Larmor radius. 

In addition to a phenomenological derivation of spectral scaling laws, summarised in 
section [5. 4[ we have given derivations of exact third-order laws for the forward-cascade 



range and for the inverse-cascade range (equations 6.11 and 6.18). These derivations are 



based upon general considerations of symmetries of the gyrokinetic equation, and the 
outcome is valid uniformly at all scales, above and below the Larmor radius. 

Combining the exact third-order results with the dual self-similarity hypothesis (equa- 
tion 8.1 ), we have been able to reproduce the phenomenological scaling laws for both the 
forward-cascade range and the inverse-cascade range (section |8]). The dual self-similarity 
hypothesis is also used to predict scaling laws for the phase-space spectrum of free energy. 
To describe the velocity-space structure, we have introduced a spectral formalism based 
on the zeroth-order Hankel transformation (section [7|. This novel spectral treatment of 
perpendicular velocity space may in general prove useful for theoretical and numerical 
applications in gyrokinetics. The works Tatsuno et a/.| ([2009, '20106 1 give strong numer- 
ical evidence of the fc-space scalings, first predicted by Schekochihin et al. (2008 2009), 



and our new predictions for the phase-space Fourier-Hankel spectrum Wg{k,p) in the 
forward-cascade range. Numerical investigation of other novel predictions of this paper, 
including the exact laws and the spectral scalings for the inverse cascade, is left for the 
future. 

9.2. Forcing and universality 

This paper explores the fundamental nonlinear processes that govern the turbulent evo- 
lution of a magnetised plasma. For the sake of both simplicity and generality, we do not 
fix the particular mechanism that drives this turbulence. It may, for instance, be driven 
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by linear instabilities that are induced by large-scale (equilibrium-scale) gradients {e.g., 
see | KronimeS|,2002J . In this case, additional (linear) terms in the gyrokinetic equation 



must be included to account for the variation in the background equilibrium; simple phe- 
nomenological arguments indicate that the additional terms do not invalidate the inertial 
range assumption at fine scales. The details of these arguments are left for future work. 

We think that the sub-Larmor scalings derived in this paper are likely to turn out to 
be universal properties of gyrokinetic turbulence. We are encouraged in this view by the 
spectra of potential fluctuations rather similar to ours having been reported by |Gorler| 
& Jenko (2008) in their three-dimensional simulations of tokamak turbulence. 



9.3. Fusion 

It is often noted that the macroscopic properties of a fusion plasma (i.e. transport of 
heat and particles) are most directly influenced by scales of the energy injection or larger 
(the Larmor scale in fusion plasmas). This has meant that inertial-range (sub-Larmor- 
scale) theory is not a focus of present-day magnetically-confined fusion research. Yet, 
understanding the phase-space cascade is central to understanding the kinetic nature 
of nonlinear interactions. Furthermore, it is a robust mechanism for dissipation in gy- 
rokinetics. Indeed, the limit of weak collisionality will be increasingly relevant in the 



scaling to larger machines. In this limit, the dissipation scale (equation 5.14) decreases, 
while the nonlinear phase-space cascade ensures that fluctuations access this scale (in 
the limit of small collisionality, the (free energy) dissipation rate approaches a non-zero 
constant, as demonstrated numerically by Tatsuno et al. 2009). The "practical" appli- 



cation of inertial-range studies may, like in fluid turbulence, be in the development of 
subgrid-scale models for large-eddy simulation and the testing of nonlinear simulation 
codes. 

9.4. Future work 

There are some issues not addressed by this work that should be given high priority. A 
more complete investigation of the gyrokinetic turbulence cascade would include more 
physical content, such as electromagnetic affects and the correct linear instability drives 
alluded to above. Electromagnetism introduces a fluctuating magnetic field into the non- 
linearity, which couples to the distribution function through Ampere's law and will en- 
rich the problem substantially, e.g., populating the sub-Larmor range with kinetic Alfven 
waves ( [Schekochihin et a/.||2009[ ). 

Effects of the third spatial dimension must be investigated. For instance, the assump- 



tion that the parallel streaming term (see equation 2.2 1 may be neglected requires careful 



justification. In the context of tokamak turbulence, a conventional assumption is that one 
may take the size of fcy to correspond to the distance along a field line from the stable 
side to the unstable side of the toroidal magnetic surface. If this assumption holds, the 
nonlinearity will dominate over the parallel streaming term for sufficiently small (per- 
pendicular) scales. In a more general context (i.e. not specific to fusion), the term is 
believed to play a special role in the Alfven wave cascade — which can be treated within 
gyrokinetic theory. In this case, it is believed that a critical balance condition is satisfied 
where the generation of finer scales in fcy ensures that the parallel term must be included 
(at dominant order) at all scales in the cascade forward ( Schekochihin et ar||2009 ). 

The linear phase-mixing process, due to the parallel streaming term, generates fine 
scales in the parallel- velocity (wy) dependence of the distribution function. The nonlinear 
simulations and theoretical work of Watanabe & Sugama ( 2004 ) demonstrate how this 



parallel phase-mixing, taken in isolation from nonlinear phase mixing, plays a role in 
the turbulent steady state. However, the question of how exactly the two phase mixing 
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processes coexist is open. It has been heuristically argued by Schekochihin et al. { 2008 1 
that the Unear phase-mixing process is much weaker than nonhnear phase-mixing at 
smah scales. If that is true, the scahng theory developed in this paper should still hold 
for the inertial range of the forward cascade in three-dimensional gyrokinetics. 

In general, the present capabilities of numerical simulations allow a fully resolved inves- 
tigation of (perpendicular) nonlinear phase-mixing or (parallel) linear phase-mixing, but 
not both. Thus, there is currently an opportunity for theoretical investigation of the full 
five-dimensional theory (that is, three position dimensions and two velocity dimensions) 
to pave the way to the full understanding of the turbulent steady state in magnetised 
plasma turbulence. 
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Appendix A. Normalisation 

The gyrokinetic equation, equation |2.2[ is cast in dimensionless variables. The normal- 
isation is summarised in the following table: 

tvth/L-)-t, x/p->-x, y/p->-y, z/L->-z, 

v/vth-)-v, ipqL/{TQp) ip, hv^^^L/{nop) ->■ h, Fov^^Jno Fq. 

We have adopted the following definitions: the equilibrium density and temperature of 
the species of interest are no and Tq; the mass and charge are m and q; the thermal 
velocity is wth — \/T\iT>Tis] the Larmor radius is p — wth/^^ where the Larmor (cyclotron) 
frequency is f2 = qBo/m. L is the reference scale length (i.e. system size), satisfying 
p/L ~ e ^ 1 for consistency with gyrokinetic ordering — no further mention of L will 
be needed. 



Appendix B. Response models for the non-kinetic species 

For the sake of generality, we treat both ion and electron (Larmor radius) scale turbu- 
lence, but with only one kinetic species. For ion-scale turbulence, we solve the gyrokinetic 
equation for ions and assume that the perturbed electron distribution function is given 
by the electron Boltzmann response (normalised to ion units): 



qToe 



(Bl) 



where e is the magnitude of the electron charge, and q is the ion charge. As an alternative 
to the Boltzmann response, we also consider the 'no response' model, Sf^ — 0. This is 
the exact two-dimensional solution, whereas the Boltzmann electron response is consid- 
ered "quasi two-dimensional" because a small but non-vanishing parallel wavenumber is 
assumed. A discussion of this can be found in Hasegawa & Mima (1978). Ion gyrokinet- 
ics with Boltzmann or no-response electrons is appropriate at all scales larger than the 
electron Larmor radius. 
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For the case of electron-scale turbulence (at the electron-Larmor-radius scale and 
smaller), the appropriate ion response model is again the Boltzmann response (nor- 
malised to electron units): 

^/.«-^^^^o. (B2) 
The choice of response model determines the parameter t in the quasi-neutrality con- 



straint, (equation 2.6): 



kinetic ions, Boltzmann electrons, 

kinetic ions, no-response electrons, (B 3) 

kinetic electrons, Boltzmann ions. 

Note that although it may be useful to think of the ion-scale turbulence {e.g., for com- 
parison with the Hasegawa-Mima turbulence), our results are more robust for electron- 
scale turbulence, as the ion Boltzmann response becomes more exact at high wavenum- 
bers whereas the Boltzmann-electron approximation will be broken for high wavenum- 
bers, approaching fcpe ~ 1. 



Appendix C. Derivation of phase-space spectra 

C.l. Derivation of equation \'8.14\ 
In the limit k p, the scaling of the spectrum Wg{k,p) can be inferred from the scaling 
of 5*2 (see equation 8.4) in the limit that i <^ iy. Combining the scaling hypothesis given 
by equation 8.1 with the solution for the scaling index, hg — 1/6 (equation 8.11a), we 
may express 5*2 as 



^2 = 14,1'/'^ 



(Cl) 



where ip is an unknown scaling function. Now, assuming that for i — and iy the 
structure function is non-zero, 5*2 (^ = 0) 7^ 0, we have ■0(0) = V'o 7^ 0. Taylor-expanding 
tjj around £ / £y = we obtain 



S2~\ev\ 



1/3 



£ 



When we substitute this into the expression for Wgi{k,p), equation 
term is negligible for fc 7^ 0. Therefore, 



8.4 



(C2) 

the leading order 



Wglik,p) 



£d£ vdv d£yJo{k£) cos{p£y 



(C3) 



Equation C3 implies the power law Wg{k,p) oc p ^/"^fc ^ . Repeating this analysis for the 



limit k <^ p completes the derivation of equation 8.14 



C.2. Derivation of equations \8 . 1 ^ and \8.1(^ 



Integrating equation 7.11 over p, we obtain (for non-zero k) 

Wgl{k)= I dpWgl{k,p) 



k 

TT- I £d£vdvJo{k£)S2{£,v,v). 



(C4) 
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Integrating equation |8.4| over k, we obtain 
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{p)= J dkWgiik,p)^^^J vdvd£^cos{pe^)S2{e = 0,iy,v). (C5) 



The scaling of 5*2 for zero gyro-centre increment i = can be obtained from equation |C 2[ 
The scahng of S2 for zero velocity increment iy is obtained analogously. These expressions 
imply equations |8.15| and |8.16[ 
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